1. IntroductionWhen a particle is dragged through a heat bath by an external force, it occasionally moves in the opposite direction to the external force. Such entropy consuming events have been shown to obey the fluctuation relation.[1,2] Fluctuation relations, which go beyond the traditional second law of thermodynamics, express the relative possibility to find entropy production and entropy consuming events in small-scale non-equilibrium systems. Its validity has been proved in a surprising range of experiments,[3–6] even extended to several non-thermal systems nominally outside their realm of applicability.[7–14]
Active matter systems consist of individual particles which can move by consuming energy and converting it into self-propulsion force.[15–20] Such a self-propulsion mechanism drives the system out of equilibrium, and is typically noisy in itself.[21] Random motions of active particles often work as “active bath”,[22,23] and provide athermal fluctuations for passive particles immersed in the system. It is an intrigue question to ask whether microscopic thermodynamic laws such as fluctuation relations hold in active bath.[24–26]
In this paper, we present a statistical analysis of angular velocity fluctuations of an asymmetric rotor driven by self-propelled robots. Our main results are as follows: (1) rotor angular velocities, from both experiments and numerical simulations, have strongly non-Gaussian distributions; (2) relative probability for the rotor to rotate in positive and opposite directions is related to the averaged angular velocity magnitude in a period of time; (3) symmetry properties in the probability density distributions are related to the detailed fluctuation relation (FR) of entropy productions at the trajectory level.
2. ExperimentOur experimental setup is shown in Fig. 1(a), consisting of an asymmetric rotor and self-propelled robots. The rotor is cut from a 2-cm-thick styrofoam sheet by a CNC foam cutter. The outer radius of the rotor is 11.5 cm, and the radius of 8 inner corners of the rotor is 8.5 cm. Part of the interior material is removed to reduce the rotor mass and finally each rotor weights around 8 g. The rotor is connected to an axis fixed to the bottom styrofoam board by two low-friction ceramic ball bearings.
Robots, serving as the heat bath in our experiments, are commercially available toys, Hexbug, whose body is 4.3 cm long and 1.2 cm wide. Each robot is driven by a vibration motor. The robot moves with intrinsic rotation and translation noises, which are approximately Gaussian.[27] As shown in Fig. 1(a), the robots are placed in a circular experiment cell, with flower-shaped acrylic boundary. Moving robots interact with the rotor through inelastic collision. Upon collision, the robots lose their velocity components normal to the rotor boundary and begin to slide along the surface. Such inelastic collisions drive the asymmetric rotor into unidirectional motion (CCW) with apparent fluctuations, see Fig. 1(b).
Motion of the rotor and robots is recorded by a high-speed camera placed vertically above the experiment arena, with a resolution of 900 pixel × 900 pixel over a field of view of 50 cm × 50cm. The frame rate used is 20 fps and each video is 1000 s long. Reflective sheets are attached to both rotor’s edge and robots to increase the contrast. We use a standard particle tracking algorithm to locate and track the reflective markers. The velocity of each marker is computed as the difference of location between two successive images. And rotor’s angular velocity in each step, Ω, is obtained by averaging the angular velocities, relative to the center axis, of the 8 markers on the rotor. We cut the raw angular velocity time series into segments of 200 s long. The mean value of the angular velocity and corresponding distribution in each time interval are stable after the beginning 200 s. We then exclude the first 200 s interval, and use remaining data for statistical analysis.
We define a time averaged angular velocity over a time interval τ,
, where ⟨ · ⟩ denotes an average over the entire time series, and P(Ωτ) denotes its probability distribution function (PDF).
A typical time series and PDF of angular velocity for a rotor driven by 15 robots are shown in Figs. 1(b) and 1(c). Both positive and negative fluctuations around the mean value are observed. In order to calculate time series of Ωτ, we divide the Ω/〈Ω〉 series into bins of length τ, and average over overlapping bins, where the center of each bin is shift from the previous one by a time difference 0.05 s, to improve statistics. The probability distribution functions of the coarse-grained angular velocity for different integration time τ = 0.1 s, 1 s, 2.5 s, 5 s are presented in Fig. 1(c). One could see that these PDFs are non-Gaussian even after long time integration. Motivated by findings in Refs. [12,14], we then examine the relative probabilities of positive and negative coarse-grained angular velocity Ωτ. We are surprised to find that the symmetry functions F(Ωτ) = log[P(+ Ωτ)/P(–Ωτ)] are almost linearly related to Ωτ for τ = 0.1 s, 1 s, 2.5 s, 5 s, with a slope kτ increasing almost linearly with τ as shown in Fig. 2, consistent with the results in Refs. [4,5,7,29,30]. The symmetry functions deviate from linear relation at large angular velocities, where the amount of data is often insufficient.
3. Numerical simulationTo investigate the necessary factors that lead to the non-Gaussian fluctuations and unusual symmetry properties, and to explore their existence in a wider parameter space, we create a simulation model composed of self-propelled rods and a rotary rotor, see Fig. 1(d). Rods are constructed by 5 overlapping spheres joined along a straight line and their sizes are chosen according to the width of Hexbug robots used in experiments. All spheres from different rods interact with each other via a Yukawa potential. The flower-shaped boundary and the rotor are constructed from a string of the same kind of spheres.
In simulations, center-of-mass r(t) and orientation angle θ(t) of a rod are controlled by the second order Langevin equations
![](cpb_29_6_060502/cpb_29_6_060502_eqn1.gif)
![](cpb_29_6_060502/cpb_29_6_060502_eqn2.gif)
where
m0 and
I0 are the mass and moment of inertia of each rod, and translational and rotational noise terms are included to capture the naturally existing fluctuations in robots’ velocity and moving directions.
[27] ξr and
ξθ are the standard Gaussian noises. The translational and rotational noises are not fundamentally coupling for dry active particles,
[28] and their magnitudes are controlled by diffusion coefficients
D0 and
DI0, which are separately determined by experiment measurements.
fex and
τex denote the external force and torque from other rods and boundary,
α0 and
αI0 are the friction coefficients. Energy is injected into the system through self-propelled force
fsp acting on the last particle of each rod and parallel to the rod body, which ensures the rod to easily slide along the rotor boundary after collision. This detail is important to repeat apparent uni-directional rotation observed in experiment. Change of rotor angle
Θ is described by the following equation:
![](cpb_29_6_060502/cpb_29_6_060502_eqn3.gif)
where
I is the moment of inertia of the rotor, and friction coefficient
αI is set to a very low value to simulate the low-friction ball bearing in experiment. The rotational noise term is excluded in Eq. (
3) because the gear is a macroscopic object and doesn’t show any spontaneous fluctuations.
Tex is the total torque due to collisions with the self-propelled rods. Parameters used in the simulation are listed in Table
1. With these parameters, our numerical model can reproduce the experimental results in angular velocity statistics (right column of Fig.
1) and symmetry functions (right column of Fig.
2). We explore the effect of the shape of the central gear, and the magnitude of self-propelling force
fsp, and find that the non-Gaussian shape and symmetry properties of the rotor angular velocity distributions are qualitatively maintained. Figure
2(e) shows that the linear dependence of slope
kτ vs.
τ is valid in simulation when
τ ≥ 0.5 s (Fig.
2(e)). Characteristic time scale of one self-propelling rod colliding and pushing the rotor is approximately 0.3–0.5 s, which is comparable to this critical time scale
τc. It requires further investigations in the future to determine whether these two timescales are fundamentally related.
Table 1.
Table 1.
![](https://cpb.iphy.ac.cn/html_resources/images/table-icon.gif) Table 1. Parameters in numerical simulations. .
I/g⋅cm2 |
1500 |
αI/N · cm · s |
0.015 |
m0/g |
7 |
fsp/N |
0.2 |
α0/N · cm−1 · s |
0.01 |
D0/N2 |
0.00002 |
I0/g⋅cm2 |
10 |
αI0/N · cm · s |
0.001 |
DI0/N2 · cm2 |
0.00001 |
| Table 1. Parameters in numerical simulations. . |
4. Entropy productionWe use the idea of entropy production[2,31] to deepen our understanding of the symmetry functions in Figs. 1 and 2. To compute entropy production along a trajectory, we need to know the transition probability between states. To this end, a discrete time Markovian model for angular velocity time series Ωτ is constructed for coarse-grain time larger than typical rotor-robot collision time (0.5 s). For example, we can coarse-grain typical data set (c.f. Fig. 1(b)) with an average time τ = 0.5 s, then calculate difference in two successive steps d Ωτ n = Ωτ(n + 1) – Ωτ n as the acceleration at time step n. After binning angular velocity time series Ωτ to 26 discretized states Ωi, we then can easily calculate PDF of acceleration for each specific state Pa(d Ωτ|Ωi). In the end, Pa(d Ωτ|Ωi) severs as a discrete time, discrete state stochastic model for rotor’s angular velocity.
Based on this stochastic description of angular velocity evolution, we can calculate trajectory-dependent total entropy production Δ Stot as defined in Refs. [2,31] for any trajectory of time length nτ. Here we directly follow the method in Refs. [32,33] since the angular velocity is of odd-parity under time reversal operation. Let x(nτ) denote a trajectory of length nτ, (Θ0,Ωτ0),(Θ1,Ωτ1),…,(Θn,Ωτ n), and x†(nτ) presents its time reversal process, (Θn,–Ωτ n), (Θn–1, –Ωτ n–1),…, (Θ0, –Ωτ0). Ωτ i denotes the coarse-grained angular velocity at time t = iτ, (n and i are integers). The probability of observing trajectory x(nτ) is given by P[x(nτ)] = p0(Θ0,Ωτ0)p(Θ1,Ωτ1|Θ0,Ωτ0) ⋯p(Θn,Ωτ n |Θn–1,Ωτ n–1). p0 denotes the initial distribution, and p(Θi,Ωi|Θj,Ωj) is a transition probability between two states, which equals to Pa(Ωi –Ωj|Ωj) given by our stochastic model. Since our system has rotationary symmetry, the position Θ0 dependence can be neglected and the above expression could be easily written in the following shorter form:
![](cpb_29_6_060502/cpb_29_6_060502_eqn4.gif)
And its reverse path
x†(
nπ) is constructed by initiating at the end of the forward process, reversing the sign of the angular momentum, and evolving backward in time. We have
![](cpb_29_6_060502/cpb_29_6_060502_eqn5.gif)
The total entropy production is defined as the logarithm of the ratio between the probabilities of the trajectory and its time-reversed counterpart Δ
Stot = log(P[
x(
nτ)]/
P[
x†(
nτ)]). The entropy production characterizes the irreversibility of the forward path. According to how the reverse trajectory is defined, we choose the initial probability of reverse path
![](cpb_29_6_060502/cpb_29_6_060502_ieqn2.gif)
as the final probability of the forward process
pn(
Ωτ n), as mentioned in Ref. [
33]. The change in trajectory entropy log(
p0(
Ωτ0)/
pn(
Ωτ n)), which equals to
![](cpb_29_6_060502/cpb_29_6_060502_ieqn3.gif)
, is naturally included in the definition. Since the rotor rotation reaches a steady state in our case, we set both initial and final distributions of
Ωτ as the steady state distribution
p0(
Ωτ) =
pn(
Ωτ) =
P(
Ωτ), which has been provided in Figs.
1(c) and
1(f).
Given a trajectory consisting of a series of n Ωτi, we first allocate each value into corresponding discrete state, and calculate the probabilities of the forward and backward processes using Eqs. (4) and (5). We can then calculate the entropy production for each trajectory. The linear fluctuation relation
is found to be satisfied for any trajectory time length (Fig. 3(a)). This result in turn proves the validity of the model we use to describe the rotor’s motion. The PDF of Δ Stot is non-Gaussian, similar to the distribution of angular velocity (Fig. 3(b)). The same technique can be applied to simulation data. Linear relation with slope 1 and non-Guassian PDF for Δ Stot are also observed in the numerical result (Figs. 3(d) and 3(e)).
5. Physical significance of stochastic entropy productionTo illustrate the significance of the computed entropy production, we now turn to the connection between Δ Stot and angular displacement ΔΘ for a specific trajectory. For a given trajectory length nτ, we calculate the average value and standard deviation of Δ Stot for all trajectories with a specific angular displacement ΔΘ. We plot Δ Stot/nτ against ΔΘ/nτ in Fig. 3(c). Here we divide both of quantities by nτ to put data from different trajectory lengths on the same plot. Note that ΔΘ/nτ is exactly the average angular velocity Ωnτ. We find that the data collapse for long trajectories (nτ⩾6 s) and an almost linear relation between Δ Stot/nτ and ΔΘ/nτ (i.e., Ωnτ) appears, with a slope factor not equal to 1. This fact shows that, in our system, the average angular velocity during a period of time can server as a direct indication of irreversibility. It also explains the FR-like linear dependency shown in Fig. 2.
Statistics of angular acceleration contains important information. As shown in Figs. 4(a) and 4(c), the average value of d Ωτ decreases almost linearly with rotation velocity Ωτ, suggesting a linear drag term in the governing equation. The variance of d Ωτ,
, also depends on the rotation velocity, as shown in Figs. 4(b) and 4(d): variance for Ωτ < 0 is apparently larger than that for Ωτ > 0. A probable reason for this dependence is a strong coupling between rotor motion and nearby robots. As shown in our previous publication,[27] rotor’s counter-clockwise rotation (in positive direction) is able to strongly regularize the positions and velocity orientations of outside robots while rotation in negative direction leads to more chaotic outside robot motion, hence larger variance in angular acceleration.
With results of the mean and variance for d Ωτ, we can write a discrete time equation for the angular velocity
![](cpb_29_6_060502/cpb_29_6_060502_eqn6.gif)
including a driving term
M, a viscous drag
γ, and a noise term with variable amplitude
f(
Ω) (
f(0) = 1), 〈
ξi ⋅
ξj〉 = 2
Dδij/Δ
t. The
M and
λ values can be obtained by fitting mean value d
Ωτ;
D can be estimated by the variance of d
Ωτ at
Ωτ = 0. By numerically solving Eq. (
6), we are able to reproduce major results presented in Figs.
2 and
3. The distribution of noise term
ξ does not have a significant effect on our results. If
f is a constant 1, equation (
6) reduces to a regular Langevin equation which has been well studied;
[14] in this case, the heat bath temperature can be found as
T =
D/
γ and the entropy production is related to the angular displacement as Δ
S =
MΔΘ/
T =
MγΔΘ/
D. Inspired by this result, we plot Δ
Stot/Δ
Θ against
Mγ/
D in Fig.
5 and observe a linear relation. This suggests that rotor dynamics in our experiments is governed by a Langevin-like mechanism with a constant drive and a viscous drag. It is possible that
D/γ gives an estimation of the effective temperature of the heat bath of active robots. Obviously, it requires more experimental and numerical data to determine the significance of such a quantify, which is beyond the scope of this paper.
6. Discussion and conclusionWe have studied angular velocity fluctuations of a rotor driven by self-propelled robots. Both experiments and simulations show that relative probabilities of observing rotation in opposite directions are linearly related to the angular velocity magnitude. To understand this observation, we constructed a stochastic description of the coarse-grained angular velocity time series and showed that the linear dependence found in angular velocity symmetry functions originates from the fluctuation relation of stochastic entropy production at the trajectory level. Both probability density functions of the entropy production and the directly measured rotor’s angular velocity are non-Gaussian, which is explained by a strong coupling between the orientation of nearby robots and the rotating state of the rotor. Our data also have shown that rotor motion bears strong similarities to Langevin dynamics with a constant driving torque and a linear drag, which explains the linear relation between the entropy production and angular displacement in a specific trajectory.